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INTRODUCTION 

Optical (van Paradijs et al. 1997) and X-ray (Costa et al. 1997) afterglows of gamma-ray bursts 
(GRBs) had influenced and still significantly affect the evolution of our understanding of the nature of 
GRBs (see e.g. reviews Postnov 1999; Meszaros 2002). 

It is now commonly accepted that the powerful energy release of the GRB itself is associated with 
gravitational energy transformation during the massive star's core collapse into a black hole or into a 
neutron star, and the broadband afterglow is believed to arise from interaction of relativistic ejecta (from 
the GRB central engine) with the surrounding medium (see e.g. Piran 2004). Also observations suggest 
that some fraction of long GRB-s is connected to unusually bright supernovae (see Woosley & Bloom 
2006 and references there). 

Optical afterglow variability at different times and on the different timescales is sometimes detected. 
In some cases (GRB 041219A, 050820A, 080319B) it was observed at the same time as the prompt 
emission, but more often it exposes itself as irregular deviations from the power-law fading of the optical 
flux at times about 10 4 -10 5 s since the GRB begins (GRB 021004, 030329 etc.). These irregularities could 
be explained by an interaction between expanding fireball and inhomogeneous outer medium (Lazzati et 
al. 2002) or by an additional energy injection (Bjornsson et al. 2004). 

Obviously, structures of circumburst matter may variate strongly , and the interaction of the 
prompt emission with these structures may cause transient phenomena lasting up to several years. 
(Bisnovatyi-Kogan & Timokhin 1997, Barkov & Bisnovatyi-Kogan 2005, Postnov et al. 2004). 

Having at our disposal a sample of 58 GRB-s with known redshifts and optical afterglows (Badjin 
et al. 2009) we had performed a search for transient features in their optical lightcurves. A notable 
part of objects was found to have irregularities, i.e. deviations from power-law fading, which looked like 



temporary rebrightenings or decrease of the fading rate. The table 1 and the figure 1 present typical 
parameters and the shape of the irregularities found. Unfortunately we were unable to obtain data of 
X-ray observations for the times when the irregularities were detected for the GRB-s listed in the table 1, 
but it is well known that X-ray afterglows usually demonstrate even more complicated temporal structure 
than optical ones (see e.g. Gehrels et al. 2009). 

Such effects (in the optical band) can be explained by the radiation of some additional energy 
either due to the "Late- Jet - Cocoon interaction" (Shen et al. 2009), or due to cooling of the matter 
previously heated by the main prompt emission. The characteristic times and variability scales in the 
former scenario seem to be significantly shorter than in the latter one. Obviously both scenarios can 
take place, but some evidences of the second one can be indicated, such as X-ray spectral lines detection 
(Postnov et al. 2004 consider GRB 011211, and Gehrels et al. 2009 also discuss some other cases) and 
presence of the "plateau" at late stages of the optical afterglow lightcurves. This article is dedicated to 
the studying of heating, cooling and radiative processes of the circumstellar matter being illuminated by 
the GRB prompt emission. 

Barkov and Bisnovatyi-Kogan in their work (Barkov & Bisnovatyi-Kogan 2005) considered a similar 
problem of single delta-like pulse Compton heating of an extended optically thin molecular cloud in terms 
of 2-dimensional hydrodynamics. Unlike their work, in our simulations we did not make any assumptions 
about the surrounding matter optical depth or about its full ionization state, we also used a continuous 
prompt emission luminosity lightcurve, took the photoionization heating into account in addition to the 
Compton scattering, and modeled the resulting thermal radiation by means of the transfer equation 
solving. Because of significant conceptual and methodological differences we recommend to read their 
paper separately, as an example of another approach, and in what follows we shall not discuss the results 
of Barkov & Bisnovatyi-Kogan 2005 or compare them with ours. 

DEFINITION OF THE PROBLEM 

As the basic concept of our model we consider a relatively dense shell around the GRB progenitor 
(the shell may be previously released by pre-GRB star some time before the explosion), which is being 
heated up to high temperatures, by the Compton scattering and by the photoionization, and then radi- 
ating its thermal energy. The goal is to calculate spectra and lightcurves of this radiation, taking into 
account non-stationary heating-cooling and ionization processes, radiative transport and hydrodynamics. 

Making the assumption that optical luminosity isotropic equivalent is to be about 10 43 erg s^" 1 
and cooling function A <~ 10 -21 erg cm 3 s _1 , one can estimate the radiating volume and obtain the limits 
on the shell geometrical parameters: V ~ Tr9 2 et R 2 h ~ 3 x lO 43 ?^ 2 ^ 4 cm 3 , where 0j et is gamma-ray 
collimation angle, R - the shell average radius, h - its thickness, ran - particle density in 10 11 cm -3 
and rj op t is the relative effectiveness of the optical emission. The 9j et angle must be of order of several 
degrees, and R <~ 10 15 — 10 16 cm, i.e. shell must be distant enough from the GRB source to let the prompt 
emission arise, but not too far to keep the photon density high enough to heat the medium significantly. 
The needed shell thickness h is therefore defined mainly by the optical efficiency. In our calculations we 
use initial h ~ 10 13 10 16 cm, n <~ 10 6 — 10 11 cm~ 3 and the density profiles both the uniform one and 
the windlike one (p oc r~ 2 ). 

The question of how such shells can be formed lies outside of scope of this work, but there is some 
evidence indicating that very massive stars can release a huge amount of mass (up to several M©) due to 
pulsational instabilities (Woosley et al. 2007). As the released shells interact with each other or with the 
outer medium, complicated structures can be formed with the different densities and the density profiles 
such as thin dense "walls" . 
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THE METHODOLOGY AND THE CODE 

BASIC PRINCIPLES 

To model the processes of interest numerically we used the STELLA multigroup radiation hy- 
drocode (Blinnikov et al. 1998) modified to reproduce the non-stationary processes of heating and the 
ionization state setting along with the hydrodynamics and the multiwavelength radiative transfer. The 
original code STELLA is described in Blinnikov et al. (1998), therefore in this section we will focus only 
on the modifications made. They are following: 

- we have added into a heat balance equation (2.11 in Blinnikov et al. 1998) a time-dependent 
term responsible for the gamma-ray heating due to the photoionization and the Compton scattering on 
free and bound electrons; 

- in those radial zones, where gamma-ray heating occurs, the state of the matter is computed by 
solving a system of non-stationary ionization balance equation, which takes into account mechanisms of 
collisional ionization, radiative ionization (the photoionization and ionization due to Compton scattering), 
radiative and diclcctronic recombination under the assumption that all ions are in ground state; in those 
zones, where heating either has not yet begun or has already ended, the ionization state is still computed 
in Boltzmann-Saha approximation (like in the original STELLA code), and a smooth connection of these 
two approximations is performed just after the heating ends; 

- special program tools have been added to take into account how the gamma-ray beaming and 
shell curvature change the lightcurves and spectra. 

One-temperature one-dimensional fluid model without magnetic field is adopted to simplify the 
calculations. It means that one has to treat the lightcurves at their early stages as averaged on the 
timescale of electron-ion energy equipartition (which in turn depends on medium parameters and model 
assumptions about electron-ion collision mechanisms, in our conditions this timescale is to be not greater 
than few hundred seconds for the Coulomb collisions only). Also in the one-temperature model we 
underestimate the early stage hard X-ray luminosity, because an amount of high energy electrons occurs 
to be underestimated significantly (following Sazonov et al. 2003, one can obtain that the electrons can 
be heated by the Compton scattering to temperatures of tens keV while in the equipartition approach 
their characteristic energy has to be three orders of magnitude lower) . 

One-dimensionality leads to underestimation of a sideways expansion of the heated matter of the 
shell into neighbouring cold regions. In our model, however, we simulate some two- and three-dimensional 
effects such as fluid velocity redistribution acceleration due to mixing (see appendix B in Blinnikov et al. 
1998) and temporal delay due to geometrical curvature. The second effect leads to temporal smoothing of 
the lightcurves and the spectral evolution, and it weakens the influence of the one-temperature approach 
on reliability of our simulations. 

These simplifications are significant, but they give us a possibility to pay more attention to the 
processes of non-stationary kinetics and radiation transfer, keeping the temporal resolution high (i.e. 
to use more timesteps without considerable increasing of the total calculation time) and using about 
hundred photon energy groups. 

INTERACTION OF GAMMA-RAYS WITH CIRCUMBURST MEDIUM 

We characterize the gamma-ray emission with the following parameters: spectral boundaries, 
a peak luminosity isotropical equivalent, a total duration, a spectral shape function and a tabulated 
lightcurve. 

The effectiveness of the radiative heating is determined in the first place by shell opacity. As the 
shell zones are being heated their state and opacity are being changed, affecting significantly the gamma- 
ray flux propagating towards the next outer zones. That is why it is necessary to keep an opacity table 
for each zone and for each gamma-ray phase moment in working memory. 
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In our calculations we determine the optical depth of the given zone as 



r z (n z , e, t p ) = h z Y, tot = h z (<r KN (e) (n fe (n z ) + n be , c (n z , e)) + Y, pl {n Zl e)) , (1) 

where n z is the zone number, h z - zone's geometrical thickness, e - gamma-ray photon energy, t p - 
gamma-ray phase time (time since the first gamma-rays crossed zone's outer boundary), cjkn{^) - Klein- 
Nishina cross-section, n/ e - number density of free electrons, nb e ,c - number density of bound electrons, 
which can be knocked out by the Compton scattering of gamma-photons (e/ s (e) > \ c , see below), E p i - 
inverse mean free path due to the photoionization process 

£> pi (n z , e) = ^2(n zl a zl (e)), (2) 

z,i 

where n zi is the number density of ions with atomic number Z and with (Z — i) electrons, a zi (e) - pho- 
toionization cross-sections for ground states of these ions (summarized over all atomic shells), computed 
using analytical fits from Verner & Yakovlev (1995) and Verner et al. (1996). 

Strictly speaking, the cross-section of the Compton scattering on bound electrons differs from 
the Klein-Nishina cross-section for photon energies less or near 1 keV (if we consider astrophysically 
typical elements), but this difference quickly disappears with the growth of frequency. That is why it is 
convenient to use one common Klein-Nishina formula (Tkn(^) for all the species instead of many different 
approximations for every ion and all its shells. The Compton scattering on bound electrons (the process 
i + "/ — > i + + j' + e~) plays an important role in our problem because its cross-section for gamma-rays is 
several orders of magnitude higher than the photoionization cross-section (i.e. the process with complete 
photon absorption: i + 7 — > i + + e~) for atomic levels with ionization potentials about few keV. 

A photon number spectral density in a given zone is calculated as 



e — h£tot(nz,e,t p ) 



n 7>e (e, h, t p ) — n 7>e (e, 0, t p )R zl (t p ) ^ ^ ) _|_ 



(3) 



here R Zt i - zone's inner radius, h - coordinate "thickness" (is counted along a radial axis from R z ,i), and 
n 7 , e (e, 0, t p ) - spectral photon density at the inner boundary of the given zone at the phase time t p . The 
parameter n Jt€ (e, 0,t p ) is determined as: 

n^,0,t P ) = L ^2^ tp) c f p) - ^Jt^ X 6XP (- E m P )Z tot (n' z ,e,t P ))) , (4) 

Z ' X Z ' X \ n' z <n z J 

where Li SO is a bolometric isotropic luminosity equivalent, F e (e) - the spectral shape function, rj{t p ) - 

gamma-ray lightcurve, F$ - normalizing factor. The sum of optical depths T' z (n' z ,e,t p ) in the exponent 

of (4) is calculated over all zones situated under the given one. 

These formulae represent the simplest approach to the gamma-ray transfer in which scattered 

photons are excluded from our consideration. This allows us to compute the gamma-ray photon field at 

each time step for about hundred energy points, without dramatic increase of the calculation time, that 

would inevitably happen if we solved the transfer equation for gamma-rays as well as for thermal photons 

(and thermal radiation Eddington factors now are recalculated only once per several tens of timcsteps). 

The assumption that the gamma-ray photons suffer not more than a single scattering allows us to limit 

the T z (n' z , e, t p ) array to t p in the phase time dimension. Thus we somewhat underestimate the energy of 

radiative heating and its duration, the more significantly the more shell's optical depth exceeds unity. 

In the approximation described above it is convenient to average the photon number spectral 
density over zone's thickness (in the following expression an explicit indication of dependence on t p and 
e is mostly omitted): 



{0)Rji h f e- fts '°' _ L lso F c (e) V (t p )k{e) R z , 2 Xtot - 2 + (2 - (g^ - fc 2 )S tot ) e -^ s '°' 

h z J {R zA + hy ~ 4tvR 2 z lC eF h z R z , 2 T, tot ' 1 ' 
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here R z .2 — R z ,i + h z - zone's outer radius. 

In these terms we can in general form determine a power of heating as: 

oo 

Q = cV z / n 7 ^(e)e^2(n a a a rj a )de, (6) 
o 

where V z is zone's volume, the sum over a denotes summarizing over all the processes of interaction of 
the gamma-rays with the matter, n a - number density of particles participating in that processes, a a 
- process cross-section, r\ a - a photon energy fraction that transforms to a knocked out electron kinetic 
energy. 

The STELLA code uses a specific heating power in units 10 12 erg s _1 g" 1 determined as 
. 10* 



/ %, e (£)EX 

Pz J 



1 - E pi (e) + a KN (e) (n fe (/.(e) - _^- s ,( e )) + n be , c {e) (/ s (e) - ^))) de, (7) 

where T e - an electron temperature (in our case it coincides with the ion one), functions f s = 
/(^f?) <r/J(e) and 9s = g(^?) CT /J( £ ) designate an averaged photon energy fraction that the pho- 
ton loses or gains in a single act of the Compton scattering, functions / and g are taken from Sazonov et 
al. (2003): 

, 3(x-3)(x + l)ln(l + 2a;) -Wx 4 + 51a; 3 + 93a: 2 + 51a; + 9 
f{x) = 8^ + 4x2(1 + 2x)3 ® 

_ 3(3a; 2 - 4a; - 13) ln(l + 2x) -216a: 6 + 476a; 5 + 2066a: 4 + 2429a: 3 + 1353a; 2 + 363a; + 39 
9(X >- 16a;3 + 8x2(1 + 2x)-> (9) 

In the equation (7) x(e) and Xc( e ) arc effective ionization potentials for the photoionization and 
the ionization due the Compton scattering respectively. The rule of their evaluation is based on the sum 
over processes (6): 

^zi^zinl (tjXzinl ^zi\zinl ^e,zinl 

x(e) = s-(i) ' ^ = • (10) 

zi 

Subscripts zinl designate the nl subshcll of the ion with atomic number z and (z — i) electrons, \ 
are ionization potentials, a - photoionization cross-sections (for ionization due to the Compton scattering 
<jkn is used for all zinl), N e ^ Z i n i - a number of electrons on the corresponding ion subshell (in the first 
equation of (10) it is implicitly contained in a z i n i), and the sum is taken over all ion subshells from which 
electrons can be knocked out by photons of energy e. 

The photoionization rate also depends on the photon number spectral density: 

oo 

P z i = c J n 7i£ (e) ((TKN(e)N e , z i tC (e) + cr zl (e)) de, (11) 
o 

where N etZitC (e) = £ N e, 

z i n i is a number of electrons in the ion zi that may be knocked out by 

n,l:Xzinl< f -fs 

photons of energy e due to the Compton scattering. 
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THE IONIZATION BALANCE AND THE STATE OF THE MATTER 

To determine the ionization state of the matter being heated, in the modified code STELLA we 
solve the following system: 



' 8n z0 

■ — 


— (Pzo + Czotie) n z o + a z on e n z i, 


for neutral atoms; 


dn zi 
~dt~ ~ 


{P z ,i-i + C z ,i-\n e ) n e n Zji _i — (P zi + (C Z i + a Zji _i)n e ) n zi + a z in e n z ,i+i, 


for partial ions; 


dn zz 

~~dT~ ~ 


(Pz, z -i + C z , z -m e ) n e n z>z -i - (P zz + (C zz + a z>z -±)n e ) n zi , 


for full ions; 


^{in zi 


) = n e , 


for electrons; 



(12) 



. zi 



where C z , are collisional ionization rates (initial state is zi), a Z i - combined radiative and dielectronic 
recombination rates (final state is zi) , z and i possess all values that are possible in a specified composition 
of the shell. Formulae, tables and subroutines, which are necessary to compute these rates, were taken 
from the web-site "Atomic Data for Astrophysics" http:/ /www.pa. uky.edu/~verner/atom. html . 

We solve the system (12) implicitly at each timestep (until the gamma-rays completely escape the 
heated zone) by replacement of the derivatives by finite differences = • The system does not 

contain hydrodynamics equations, so all n Z i have to be scaled with the total density variation in the zone, 
which is determined separately. But system's solution affects the heating power (7) and thus provides 
some kind of feedback on medium hydrodynamics. 

As an iteration process convergence criterion we use zone's gamma-ray opacity (for several energy 
groups), namely the value (1 — e~' lzStot )/£ tot , because during the heating phase it is of our primary in- 
terest. If the opacity is changed significantly during the iterative solving of (12), then the photoionization 
rates (11) are to be recalculated. 

Using the system's solution the photon density (of the gamma-ray coming to the neighboring 
zones), the internal energy density of particles, the pressure and their partial derivatives over the mass 
density and the temperature are calculated. 

Before the heating begins in the zone and some relaxation time (set manually, now it is 10 seconds) 
after it ends, the state is calculated in the Boltzmann-Saha approximation, like in the original STELLA 
code. During the relaxation time a smooth connection of the two solutions (the Boltzmann-Saha solution 
and the system (12), but with P Z i excluded) is performed. 

INTERACTION OF RELATIVISTIC EJECTA WITH THE CIRCUMBURST MEDIUM 

When studying processes which occur near the GRB central engine it is also important to pay 
attention to how the circumburst medium is affected by the relativistic ejecta that generates the gamma- 
rays. If a matter structure is situated at the distance R from the gamma-ray emission source, then the 
ejecta will reach it At « i?/(2cF 2 ) seconds later after the gamma-ray forward front (here L stands for 
an ejecta Lorentz-factor). Given R w 10 16 cm and r w 20 -j- 30 (see Zhang & MacFadyen 2009), the 
temporal delay is to be about At rs 180 -=- 400 seconds. 

On the one hand, the non-relativistic STELLA code does not allow to model such essentially 
relativistic processes properly, but on the other hand the short delay At does not allow to neglect them. 
We can partially avoid these troubles without complete rewriting of the code by calculating in the non- 
relativistic formalism how the matter is penetrated by "quasi-ejecta" i.e. an additional shell moving with 
the speed of light at the distance cAt behind the gamma-ray forward front and having the same isotropic 
kinetic energy (M q c 2 /2) as the realistic relativistic ejecta (« Fmc 2 ). In this case the "quasi-ejecta" mass 
M q plays role of the Lorentz-factor. 

Clearly, this replacement is quite crude, that's why we use the "quasi-ejecta" approach only to 
determine a qualitative difference of properties of the matter heated by the radiation and the matter 
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heated both by the radiation and the kinetic shock. In discussion we also describe some features of such 
modeling of the relativistic ejecta by means of the non-relativistic code. 



OUTPUT PARAMETERS AND RESULT ANALYSIS 

The main studied parameters are luminosities and fluxes in different spectral bands (the bolometric 
one, U, B, V, R, I, SXR for Soft X-Rays 0.1-2 kcV and XR for X-Rays 2-10 kcV). A thermal radiation 
luminosity spectral density is also modeled. 

For a terrestrial observer the flux is equal to 

^ Liso 



^Dl hot 



(13) 



here L iso is an isotropic luminosity equivalent, and D phot is a photometric distance. The term isotropic 
equivalent designates such a luminosity which the object would have if it radiated in all directions the 
same energy per solid angle as in the direction of the observer's detector (Ld), so 

here fi^ is detector's solid angle as it is seen from the source. This is the estimation of L iso that we can 
obtain from our observations. Thus we need to define its calculational analog to be able to compare our 
simulation results with the observed parameters correctly. 

To compute L iso we have to know a radiation pattern of the source and a brightness distribution 
over its surface. There are several additional features that complicate the problem: the gamma-ray 
collimation (which cuts-off a relatively small area of the spherical shell), possible noncoincidence of a 
jet axis with a line of sight, the temporal delay of the arrival time of the photons coming from different 
points of the shell caused by shell's curvature (St oc (1 — cos(9))R/c) . 

If the jet half-opening angle is 9j et , the off-axis angle is 8 oa (see fig. 2), and spherical coordinates 
of the shell points relative to its centre are parametrized as (8, ip), then 

L d (t) = J \(t'(8,t)) J g(8',ip')<m'dS, (15) 

S fid 

here (8' , ip') angular coordinates relative to a surface element dS centre and to its normal, X(t'(8,t)) = 
Ltot(t'(8,t))/4:TrR 2 (t'(8,t)) is a luminosity surface distribution at a time t' = t — 6t(6,t), L tot - the full 
shell luminosity at its outer edge (R), a coefficient g(8' , ip') oc cos(8') sets an angular distribution of the 
surface clement (dS) radiation, integral over dS is taken on the shell outer edge, and integral over dU,' is 

2W2 

taken over detector's solid angle. A normalizing condition for g is J J g(6' , ip 1 ) sm(6')d6'dtp' = 1, and 

o o 

it gives g(8',ip') = cos(8')/ir. Then, taking into account an infinitcsimality of the solid angle ild (i.e. 
g(8',ip') are constant inside it), the integral over dfl' can be replaced by cos(8)Qd/^- And taking into 
account that dS = R 2 (t'(8,t)sm(d)ddd<p, one can obtain 

L iso (t) = 4tt^ = - [ L tot (t'(8,t))(M8)-M8))cos(9)sm(9)d8, (16) 



here (<pi(8) — <Po(@)) indicates what part (in units of angle) of an elementary ring sm(8)d8 (dashed 
concentric rings on the left panel of the figure 2) is situated inside a gamma-ray heated spot of the shell, 
besides 8 this difference depends also on 8j e t and 8 oa . 

If the thermal radiation photon emitted from the polar angle 8 is detected simultaneously with the 
"on-axis" photon emitted at time t, then emission time t'(6,t) (by observer's watch) of the former one 
can be found from the equation 

, R(t')(l -cos(0)) 

t' + — — y —^=t. (17) 
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It is clearly seen that cos(#) sm(8)d8 = —d(cos 2 (8))/2 and therefore it is convenient to choose 
partition over 8 providing a uniform grid of value x = cos 2 (#) points. 

The most principal output of photometric calculations with the STELLA code is a luminosity 
spectral density table L to t,v for different time (up to several thousands moments) and frequency points. 
Calculations using formulae (16)-(17) are carried out for each frequency point to obtain an array of 
Li so ,v{t)- The luminosity lightcurves in different filters are computed by intergating of L iso ^{t) multiplied 
by different transmission curves over different spectral ranges. 

Also the "red shifted" array L iso ^(t) (time values are also changed in it) is used in calculations 
of fluxes and stellar magnitudes. A cosmological model with H = 73.5 km s _1 Mpc -1 , = 0.76, 
fip = 0.24 is used and no supposition about an extinction either in a host galaxy or in the Milky Way is 
made. 

NUMERICAL MODELING 

Using the methodology described above we carried out a set of calculations for several dozens of 
shells, which differed in geometrical configuration, density and temperature. Parameters of the most 
typical of them are presented in a Table 2. 

The incoming gamma-ray radiation was set by three FRED-pulses (Fast Rise - Exponential Decay) 
with a characteristic fading time of 1.3 seconds, a peak isotropic luminosity equivalent 3xl0 53 erg s~ x 
(see fig. 3) and a Band spectral shape function (Band ct al. 1993) 



with eo = 300 keV, a — 0.9 and f3 = 2.001 in an energy range 1 keV - 30 MeV. Angles were taken 
6 j e t = 10°, 8 oa = 3°. The peak luminosity was taken a bit higher than that of the real GRB-s with 
the afterglow irregularities (see the table 1), because it represented the radiation that had not yet been 
absorbed, therefore it was reasonable to add some luminosity reserve (but not great, so our accepted 
luminosity was still of typical value for GRBs). 

In calculations with the "quasi-ejecta" model, we took the mass of the "quasi-ejecta" M q = 0.5M Q 
which corresponded to an isotropic kinetic energy Eki n ,iso = M q c 2 /2 = 4.5 x 10 53 erg, i.e. to the value 
that relativistic ejecta should possess at a distance of ~ 10 16 cm from the central progenitor. Initially 
the "quasi-ejecta" shell was placed 200 light-seconds behind the gamma-ray radiation front. 

Results of calculations for typical models (see table 2) are presented in figures 4-16. 

DISCUSSION 

THE RADIATION-TO-MATTER INTERACTION 

First of all it should be mentioned that the use of the non-stationary system (12) for the ionization 
balance calculation leads to a remarkable difference from either fully ionized matter approximation or 
Boltzmann-Saha solution. The figure 4 shows the dependencies of the optical depth (in gamma-rays) on 
the photon energy calculated in different approximations. It is seen (fig. 4a), that when the system (12) 
is used, the matter quickly (much less than 1 second) becomes fully ionized, but at the very early stages, 
however, it has much higher optical depth (and therefore absorbs much more energy). In the Boltzmann- 
Saha approximation the matter temperature occurs to be not sufficient to ensure the full ionization and 
thus the shell optical depth and the absorbed energy are significantly overestimated. 

In the figure 5 spectra of the incoming and outcoming gamma-radiation in its different lightcurve 
phases. A considerable deficiency of low-energy photons at the early stages of the prompt emission is 




(18) 
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predicted (but due to the cosmological redshift this feature should fall into X-Ray band and probably 
should not be observed by gamma-ray detectors). It is also seen that outcoming gamma-ray spectra 
become harder than incoming ones and that luminosities may significantly decrease. 

The figure 6 shows temperature radial profiles at different moments during the shell heating for 
the "wall" model with 100 radial zones, i.e. temporal resolution is w 17 s. Sharp peaks seen in profiles 
for t = 167s and t — 1670s correspond to those zones where the heating has just begun. 

Geometrical extents and column density of shells are also of great meaning. A combination of these 
two parameters determines shapes and amplitudes of the lightcurves in a larger measure than the mass 
density or the total mass of the shell. Depending on the geometrical thickness and the matter column 
density there can arise shock waves in the shell, coming out from the inner layers that have absorbed 
more energy. In the "ball" model the windlike density profile the inner zones are the most opaque (in 
spite of that they have smaller mass than inner zones of other models) and they intercept the main 
part of gamma-radiation, so the Shockwave arising there heats the outer zones more effectively than the 
radiation itself. That is why the main maximum of a bolometric luminosity lightcurve takes place not 
just after the heating ends but near the moment of the shock breakout (see fig.7d). 

If the inner layer absorption is not so high or the shell has not large thickness (e.g. like "the 
wall" does), then the matter is heated mainly by the gamma-radiation. In such a case the total radiative 
cooling time appears to be shorter, the peak luminosity is higher and more energy is radiated in X-Rays 
for the first few hours (e.g. see fig. 7a). In intermediate cases both the cooling after the radiative heating 
and the shock breakout are clearly seen. 

A temporal behavior of the thermal radiation is also affected by the shell curvature. This leads 
to a smoothing of the lightcurve on the timescale of 5R/c, where 5R is a typical variation of distance of 
an emitting point from the image plane. Actually, "the wall" has its total thickness of 5 x 10 13 cm, but 
because of the shell curvature the delay of a large angle radiation is about St = R(l — cos 0j et )/c « 5000 
seconds. This effect is clearly seen in the X-Ray lightcurves (figure 8a, figure 9). 

OBSERVATIONAL DETECTION POSSIBILITIES OF THE MODELED EFFECTS 

Along with a development of a computational tool for numerical modeling, it is important to 
understand how results of its predictions can be detected observationally. The studied thermal effects are 
more likely to expose themselves as lightcurve irregularities against the power law synchrotron afterglow 
background. To reveal in what models, at which times and in what spectral bands the thermal effects 
can be observed, we draw in our graphs characteristic regions containing the main part of observed real 
afterglow points. To plot R-filter afterglow regions (both for luminosities and for fluxes corrected for 
the Galactic extinction) we took data from a database of the optical afterglow observations used in work 
Badjin et al. 2009. X-ray afterglow regions in a range 1-30 keV (also corrected for the Galactic extinction) 
were defined based on review of Gehrels et al. 2009. 

It follows from a comparison of thermal X-Ray lightcurves (figure 8) and characteristic afterglow 
regions that the thermal effects in X-rays can show themselves as either rcbrightcning (fig. 8a) or 
considerable deceleration of the total flux fading - both followed by an abrupt growth of the fading 
rate and then return to the power law (fig. 8b). Characteristic times, when these features appear, are 
about tens minutes after the GRB prompt emission begins, and characteristic durations are of few hours. 
Obviously, the dimmer the power law afterglow is, the clearer the thermal effects can be seen. Their 
luminosity and relative visibility also depends on the shell properties (see the fig. 8 and the table 2). The 
column density and the matter distribution must provide a relatively uniform profile of the heating power 
(i.e. not to capture the main part of radiation near the shell inner edge), but on the other hand an enough 
amount of heated matter must be seen by the observer. Probably models like "the wall" are the closest 
ones to a satisfaction of these opposing requirements. 
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In the figure 10 there are shown the calculated lightcurves of near-earth fluxes in ranges 0.1-2 
and 2-10 keV. Here no extinction in the Milky Way or in a host galaxy was taken into account. It 
makes difficult to compare these lightcurves with observed ones properly, so we did not show afterglow 
characteristic regions for the X-ray fluxes. However it is still possible to obtain ideas about unabsorbed 
fluxes and their temporal behavior. 

The R-band characteristic afterglow regions are drawn in the figures 7 (for luminosities in a proper 
system), 11 and 12 (for fluxes observed from earth). The optical luminosity lightcurves suggest that the 
thermal radiation is probably too weak and can be detected only at very late stages when a non-thermal 
afterglow radiation fades away. But in fluxes the modeled effects seem brighter and expose themselves 
considerably earlier (cf. the fig. 7a and the fig. 11a). This is caused by the fact that the synchrotron 
luminosity spectral density decreases with frequency (near the 10 15 Hz), while the modeled thermal 
luminosity spectral density grows, thus, after the cosmological redshift, in the same observer's band falls 
dimmer parts of the synchrotron spectrum and brighter ones of the thermal. This may compensate 
or even overcome the common decrease of flux with growth of photometric distance. That is why the 
probability of the optical detection of the thermal "irregularities" increases significantly with redshift. 
However, the most common way for the thermal radiation to become apparent is the decreasing of the 
total afterglow fading rate or "plateau" at late afterglow stages (after few days since GRB detection), 
and at the latest stages even "bumps" may be produced, if the thermal radiation is bright enough. 

It is worth to be mentioned that because of the redshift the luminosity estimations presented in 
the table 1 do not correspond to a proper system R-band (but only to that spectral ranges falling into 
the terrestrial R-band), and it is not completely correct to compare them with the calculated R-band 
luminosities directly. 

We can also suppose that we should look for observational detections of the thermal effects amongst 
GRB-s with high redshifts. And it is interesting to note that GRB 050904 (z=6.29) optical afterglow 
lightcurve do show a certain deviation from the power law resembling our calculated flux lightcurve (for 
z=6.29) in the "wall" model. Similar feature is seen in the lightcurve of the GRB 090423 afterglow 
(z=8.2). 

Besides relatively dense shells the case of a significantly more rarefied cloud was studied (the 
"cloud" model in the table 2, see also figures 14 and 15). Although its mass is several tens times lower 
than that of the other shells, its optical luminosities and fluxes are nearly of the same order. I.e. if the 
cloud had more radiating volume (with the same density), this would increase its calculational optical 
luminosity to a value about that in the table 1. The cooling time would be also increased. But such a 
calculation requires much more radial zones in the grid , and therefore much more computational time, 
so we have not yet performed it. 

THE INFLUENCE OF THE RELATIVISTIC EJECTA. PROSPECTS OF A FURTHER CALCULA- 
TION DEVELOPMENT. 

The preceding discussion was concerned with the effects that would take place, if the circumstellar 
matter structures were heated by the gamma-ray radiation only. This is satisfied either during the first 
10 2 — 10 3 seconds (for the observer on earth ) since the GRB begins, or in the case when the shell is very 
large and the relativistic ejecta slow down in its interior taking no significant effect on the main bulk 
of matter. In other cases the radiatively heated matter will be affected by the ejecta and relativistic or 
subrelativistic motion will arise in it. 

We tried to estimate the influence of the relativistic ejecta on the thermal radiation processes, 
as far as it was possible to do this with the non-relativistic code. It had occurred that even quite 
simplified calculations in the "quasi-ejecta" approach took much time: the high temporal resolution (less 
than 10~ 4 s) was needed during a long time (more than 10 4 s) because of high relative velocities of 
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matter elements. We avoided this trouble partially by artificially increasing a mixing parameter (see 
an appendix B in Blinnikov et al. 1998), to increase an effectiveness of a velocity redistribution. But 
this is quite speculative, because the realistic mixing parameter values are to be determined from 3D- 
calculations which yet were not carried out. Nevertheless, we believe, we can do some qualitatively correct 
conclusions. 

The figure 16 shows for a purpose of comparison the luminosity lightcurves in the "bolometric" 
(1-50000 A), X-ray and optical (the U-filter is taken because it will fall into red and infrared regions 
due to the redshift) spectral regions for the "wall" model and for three different types of heating: the 
radiation only, the "quasi-ejecta" only and the radiation combined with the "quasi-ejecta" (but in the 
last two cases the shell temperature of 3000 K was taken, to let the lightcurves to reflect the heating 
dynamics more appreciably). As expected, the radiative heating dominates during early stages, and the 
luminosity due to a kinetic heating becomes equal to the luminosity due to radiative heating nearly after 
1500 seconds (for a near outer observer unaffected by redshift) since the prompt emission begins. After 
that the shock heating luminosity exceeds the radiative heating luminosity, but not more than three 
times, and then they are equal to each other again until near ~ 10 seconds the radiative heating energy 
reserve is mainly exhausted, and the radiative heating bolometric luminosity tends to fade quickly. Soft 
X-ray lightcurves demonstrate the similar behavior. 

In the 2-10 keV range the emission caused by the radiative heating exceeds the kinetic heating 
emission considerably (more than 6 times) until it begins to fade abruptly. On the contrary, in the U-filter 
kinetic energy radiation dominates nearly all the time (in the fig. 16 one should see the difference only 
between the "kinetic" and "combined" optical lightcurves, because the "radiative" one is calculated for 
the initial temperature of 12 000 K, and therefore its initial level is nearly (Xi/T 2 ) 4 = 256 times higher, 
as it is clearly seen in the fig. 16b). 

One can conclude that the thermal energy of the gamma-ray conversion is radiated mainly in higher 
frequencies and for shorter time, than that gained by the interaction with the "quasi-ejecta" . The relative 
increase of the radiation duration due to the "quasi-ejecta" influence is to be especially pronounced for 
geometrically thin shells (the duration increases by an order of magnitude), while for more extensive 
shells the difference is not so significant (we also carried out the calculations with the "quasi-ejecta" and 
the radiative heating for the "thin layer" model). If the shell is large enough, then the "quasi-ejecta" will 
be quickly decelerated in its interior and lag behind the gamma-ray radiation front considerably, and the 
bolometric lightcurve will exhibit a strong secondary maximum at a few days after burst's beginning. The 
luminosity spectral density maximum will lay in a soft X-ray or far ultraviolet region at that moment. 

Thus, generalizing, we can say that, though the "quasi-ejecta" affects the lightcurves significantly, 
it is still necessary to take into account the effects of the gamma-ray-to-heat conversion. It is also obvious 
that the relativistic ejecta colliding with a dense structure (several M ) should become non-relativistic 
rapidly. If so, the replacement of its Lorentz-factor by a "quasi-ejecta" effective mass will lose its original 
meaning of a kinetic energy storage, because the Lorentz-factor will to drop to T « 1, but the "quasi- 
ejecta" mass will hold unchanged. While the low rest mass relativistic ejecta will dissipate its kinetic 
energy in the first zones and will turn into a relatively weak non-relativistic disturbance (itself, but not the 
shocked matter), the heavy "quasi-ejecta" (M q — 0.5M Q ) will still influence the matter significantly. I.e. 
they strongly differ in a radial profile of an energy exchange effectiveness with the medium. That is why 
there are some reasons to believe, that the "quasi-ejecta" approach changes the lightcurves excessively, 
with respect to what could be expected from the real relativistic ejecta, at least in time domain. 

Certainly, to describe correctly the processes taking place when the circumstcllar matter is affected 
by the radiation and the relativistic ejecta of the GRB, one needs for calculations in the framework of 
relativistic hydrodynamics and radiative transfer. This sets quite a number of challenges. We will now 
point briefly several of them. Firstly, the nature of the ejecta itself is still unclear, e.g. whether it is 



11 



dominated either by baryonic load or by an electro-magnetic field. It is crucial for our understanding 
how the ejecta interacts with the medium. Secondly, a more sophisticated model of medium structures 
is required. There obviously must be some amount of matter between the central object (the GRB 
progenitor) and the shell, and this probably should change the outflow dynamics somehow. Also, the 
shells are likely to have (and, probably, they do) a complicated 3-d structure: e.g. "windows" (regions of 
rarefaction) that may let the gamma-rays to go through relatively unabsorbed, or higher density clumps 
that intercept an excessive portion of energy of the radiation and the outflow. Thirdly, it is necessary 
to use more precise approximations for the problems of gamma-ray transfer and particle kinetics (for 
example it is useful to solve a kinetic equation at least for electrons during the first seconds of heating in 
the phase time scale). Finally, because the afterglow emission is essentially non-thermal, it is desirable to 
have possibility to calculate magnetic fields for the purpose of self-consistent modeling of the synchrotron 
emission and the magnetic field influence on kinetics of particles (e.g. on an energy redistribution time). 
We believe that the experience described in this work will be useful in a more sophisticated calculational 
project dedicated to the phenomena taking place in vicinity of the GRBs. 

CONCLUSIONS 

The main results of our work are following. Using the radiation-hydrocode STELLA as the basis, 
the computational tool was developed, that allows, under reasonable assumptions, to model the processes 
of time-dependent heating and matter state evolution along with the hydrodynamical and radiative 
transfer processes. 

As one of possible applications the problem of radiation of the dense circumstellar medium struc- 
tures being heated by the GRB prompt emission is considered. The spectra and lightcurves obtained 
from the calculations allow us to assert that modeled processes can be responsible, at least in some cases, 
for the irregularities and plateaus, that are seen in several GRB optical and X-ray afterglow lightcurves. 

The modeling of the relativistic ejecta influence on the features of thermal radiation of the matter, 
previously heated by gamma-rays, reveals that it is important to pay attention to both these kinds of 
heating and also sets the problems that must be solved to do this correctly. 
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TABLES 



Table 1. Typical parameters of optical lightcurve deviations from the power law fading. 



GRB 


z 


logL^ 
erg s _1 


R,peak 

erg s _1 


rpC 

peak 

days 


J 90 

days 


^burap 

erg 


000301C 


2.0335 


52.39±0.14 [15 _ 150] 


1.45 x 10 44 


1.276 


0.9 


(6.68 ± 1.7) x 10 48 


020124 


3.198 


52.72±0.03 


3.73 x 10 43 


0.356 


0.215 


(1.9 ±0.7) x 10 50 


021004 


2.3351 




5.73 x 10 44 


0.024 


0.15 


(1.2 ±0.5) x 10 49 


030328 


1.52 


52.33±0.07 


1.28 x 10 44 


0.115 


4.09 ? 


(3.1 ±1.5) x 10 48 


030429X 


2.65 


53.05±0.07 


2.89 x 10 43 


0.751 


0.68 


(6.7 ± 4) x 10 48 


060206 


4.048 


52.49±0.03 


4.33 x 10 43 


0.111 


0.025 


(6.2 ±2.2) x 10 48 



a Maximal luminosity isotropic equivalent in a photon energy range 1 keV - 10 MeV, but for GRB 
000301C (for which only a single power law spectrum was reported) L 1 corresponds to 15 - 150 
keV 



b Maximal isotropic optical luminosity deviation from the power law fading in a proper system 
frequency range corresponding to the earth R-filter. 

c Time since GRB beginning (in a proper system timescale) when the maximal deviation occurs. 
^ Duration of the deviation while 90% of its energy is being emitted 
e Total optical energy of the deviation. 



Table 2. Initial parameters of the most typical shell models. 



Reference 
Designation 


Rin Rout 

10 16 cm 


M b 

M 


n bar 

10 10 cm~ 3 


N d , 

col 

10 24 cm- 2 


10 3 K 


Density 
profile 


Composition 


"wall" 


1-1.005 


5 


9.5 


4.75 


12 


uniform 


WBH * 


"thin layer" 


1-1.2 


10 


3.32-4.76 


7.93 


3 


windlikc 


WBH 


"thick layer" 


0.5-1.1 


10 


0.237 


14 


3 


uniform 


WBH 


"ball" 


0.1-1.1 


2 


0.016-1.84 


16.7 


10 


windlikc 


WBH 


"cloud" 


0.1-2 


0.1 


0.000357 


0.0687 


10 


uniform 


solar 



a Outer and inner radii of the shell 



Total mass of the shell 
c Baryonic number density 
^ Baryonic column density 

e Initial temperature of uniform shells or interior layers of windlike shells 

f Composition is taken from Woosley et al. (2007) (WBH stands for Woosley, Blinnikov, Heger) 



14 



ILLUSTRATIONS 
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Figure 1: The typical appearance of the afterglow lightcurve "irregularity" (GRB021004, Holland 
et al. (2003)). 
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Figure 2: Schematic pictures of the emitting region as it is seen "from observer's point" (a) and 
in section(b), bold line depicts shell's outer surface. The GRB progenitor is in point S. SA is a jet 
axis, SO — line of sight. 
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Figure 3: Lightcurve of the incoming prompt emission. 
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Figure 4: Optical depths of shells at different phase time moments. Panel (a) — the "wall" model 
with the system (12) used to calculate the state of matter; panel (b) — one of early calculations for 
the "thin layer" with Boltzmann-Saha approximation used, where photoionization and Compton 
scattering affect only the heating power. It is clearly seen that in the latter case the temperature 
is insufficient to fully ionize heavy elements, but when the photoionization is "turned on", almost 
all electrons instantly become "blown away" from their atoms by the radiation. And only at the 
last stages, when the gamma-ray photon density fades considerably, a partial recombination takes 
place. 



16 





\ 


(a) 


— ■ — initial 

— • — sec 

- » - 1 sec 

13 sec 

15 sec 


H rf 

! 7 

f 1 

I 1 








♦ J 

; / I 









e, [keVJ 




Figure 5: Spectral shape of the gamma-ray emission coming into (initial) and passed through the 
shell at different phase times for models "wall" (a) and "thick layer" (b). 




Figure 6: Temperature profiles in the "wall" model (provides the best spatial resolution, 16.67 
light seconds) taken in different times. 
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Figure 7: Optical (UBVRI), X-ray (SXR for 0.1-2 keV, XR for 2-10 keV) and "bolometric" (1- 
50000 A) luminosity lightcurves for the following shell models: "wall" (a), "thin layer" (b), "thick 
layer" (c), "ball"(d). The time is counted with respect to the prompt emission start. The subtle 
dotted lines represent the synchrotron power law afterglow characteristic region in the proper 
system R-filter. The dashed lines depict the sum of the lowest edge of the region and the thermal 
emission in R-band. 
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Figure 8: X-ray (2-10 keV, XR), soft X-ray (0.1-2 keV, SXR) and "bolometric" luminosity 
lightcurves for the following shell models: "wall" (a), "thin layer" (b), "thick layer" (c), "ball"(d). 
The subtle dotted lines represent the synchrotron power law afterglow characteristic region in the 
proper system 1-30 keV energy range. The dashed lines depict the sums of the region edges and 
the 0.1-10 keV thermal emission. 
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Figure 9: X-ray 2-10 keV flux lightcurves as they seen from different redshifts. No extinction taken 
into account. 



20 




Figure 10: Soft X-ray 0.1-2 kev flux lightcurves for the following shell models: "wall" (a), "thin 
layer" (b), "thick layer" (c), "ball"(d), as they seen from different redshifts. No extinction taken 
into account. 
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Figure 11: Optical (UBVRI) lightcurves for the following shell models: "wall" (a), "thin layer" (b), 
"thick layer" (c), "ball"(d), as they seen from redshift z=l.l. The subtle dotted lines represent the 
synchrotron power law afterglow characteristic region in observer's R-filter (i.e. it is somewhat 
different from those in GRB's proper system). The dashed lines depict the sum of the region edges 
and the thermal emission in R-band. 
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Figure 12: R-magnitude lightcurves for the following shell models: "wall" (a), "thin layer" (b), 
"thick layer" (c), "ball"(d), as they seen from different redshifts. The subtle dotted lines represent 
the synchrotron power law afterglow characteristic region in observer's R-filter. 



23 




24 



1 1 1 1 1 1 1 1 1 1 1 


— ■ — Osec 




—•—1 sec 




15 sec 






■ ■ ■ ■ ... 1 


\(a): 



IOC 1000 

e, [keV] 





e , [keV] 



(d) 



-Bol 
SXR 

-XR 



time, days 



time, days 



Figure 14: Optical depths (a), gamma-ray spectra (b), luminosity lightcurves (c), X-ray luminosity 
lightcurves (d) for the "cloud" model. 
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Figure 15: Optical magnitude (a, b) and soft x-ray flux lightcurves, as they seen from different 
redshifts, and thermal emission spectra (d) for the "cloud" model. 
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